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Towards an exact treatment of exchange and correlation in materials: 
Application to the "CO adsorption puzzle" and other systems 
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It is shown that the errors of present-day exchange-correlation (xc) functionals are rather short 
ranged. For extended systems the correction can therefore be evaluated by analyzing properly chosen 
clusters and employing highest-quality quantum chemistry methods. The xc correction rapidly 
approaches a universal dependence with cluster size. The method is applicable to bulk systems 
as well as to defects in the bulk and at surfaces. It is demonstrated here for CO adsorption at 
transition-metal surfaces, where present-day xc functionals dramatically fail to predict the correct 
adsorption site, and for the crystal bulk cohesive energy. 
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Electronic structure theory is the base for a multiscale 
modeling of materials properties and functions (see e.g. 
Ref. [1]). Obviously, if the needed accuracy is lack- 
ing at this base, there is little hope that accurate pre- 
dictions can be made at any level of modeling that fol- 
lows. For polyatomic systems, density-functional theory 
(DFT) with present-day exchange-correlation (xc) func- 
tionals has proven to be an excellent technique for cal- 
culations at this electronic-structure base. However, it is 
not as good for certain types of binding interactions. Ac- 
curate treatments of strong electronic correlations, van 
der Waals interactions, and molecular dynamics for elec- 
tronically excited states represent unsolved challenges. 
Besides numerical approximations (e.g. the basis set, 
possible use of the pseudopotential approximation, etc.), 
that good theoretical work is typically scrutinizing, a sat- 
isfying test of the quality of the xc functional was not 
possible for bigger systems, so far. Sometimes the results 
obtained with different functionals have been compared, 
and when they agreed this was taken as indicator for reli- 
ability. Though this was the best possible approach, it is 
neither safe nor justified. Exchange-correlation function- 
als are typically built on the homogeneous electron gas 
[the local-density approximation (LDA) Q], adding cor- 
rections while ensuring consistency with known sum rules 
[e.g. the generalized gradient approximation (GGA) Q], 
or they are constructed to reproduce certain data of some 
small molecules (e.g. the B3LYP functional [4|]). There 
is no systematic expansion in terms of successively de- 
creasing errors, and, there is no proof that, e.g. the GGA 
will always work more trustfully than the LDA. On the 
other hand, for wavefunction methods, several promising 
concepts exist for better xc treatments also for extended 
systems (see Ref. p, Iff and references therein) . However, 
these are not yet efficient, in particular for metals. 

It is often argued that the xc error is largely canceled 
when total-energy differences are studied, and that the xc 
approximation affects the geometry only by little. How- 



ever, this is generally not correct. A key example for the 
xc problem is the low-coverage adsorption of CO at the 
(111) surface of Pt or Cu, where the LDA as well as the 
GGA predict that the molecule adsorbs in the threefold- 
coordinated hollow site. Experiments, on the other hand, 
show undoubtedly that the adsorption site is in the one- 
fold coordinated top site. Obviously the conclusion based 
on DFT-LDA/GGA is even qualitatively incorrect; and 
when comparing the calculated energies of the two sites 
in question, the error in the energy difference is indeed 
significant: In the LDA it must be larger than 0.4 eV. 

The comprehensive study of CO at Pt(lll) by Feibel- 
man et al. [7[ set the ball rolling, showing that when 
properly realizing all technical aspects of the calculations 
(e.g. basis set, supercell, cluster geometry) the LDA and 
GGA put the molecule at the wrong adsorption site. For 
several other close packed transition metal surfaces, in- 
cluding Cu(lll), the situation is analogous (e.g. [8j). 
Shortly after the paper by Feibelman et al. [?J it was real- 
ized that the wrong site preference of CO may be related 
to the fact that DFT-LDA and GGA inaccurately de- 
scribe the CO-molecule's chemical bond. This was then 
expressed in terms of the HOMO-LUMO splitting (the 
5cr and 2ir* energy levels) or correspondingly the singlet- 
triplet excitation, (jl. HoL 11, 12, 13] An upshift of the one 
electron 2tt* level makes a partial charge transfer from 
the 5<T to the 2ir* orbital more difficult and therefore 
stabilizes the CO bond. Indeed, such (semi-empirical) 
upshift brought the CO molecule to the experimentally 
known site. Thus, the problem appears to be understood 
(to some degree), but altogether the situation is unsatis- 
factory. A first-principles theory should provide a reliable 
answer, and an add-on, that is triggered by a disagree- 
ment with experiment, questions the usefulness of the 
whole self-consistent procedure. We also note that the 
full correction of the LDA/GGA functional will not just 
shift the CO 2ir* orbital to higher energies; it will also 
modify the substrate d-states and thereby the substrate- 
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adsorbate bonding in various terms. Below we will take 
the "CO adsorption puzzle" as our main example for a 
systematic, non-empirical approach to correct xc errors. 
And we will also discuss a fundamentally different case, 
namely the bulk cohesive energy. 

For small systems (say up to 20-50 atoms) high-level 
quantum chemistry methods or the quantum Monte 
Carlo approach can be employed to obtain accurate re- 
sults for the exchange-correlation energy. We will show 
that calculations on such small clusters are sufficient to 
evaluate the DFT-LDA error of extended systems. We 
will concentrate on CO/Cu(lll), CO/Ag(lll), and Cu 
bulk as for these systems relativistic effects are small. 
We use accurate full-potential augmented plane wave 
(LAPW/APW+lo) DFT calculations Q with LDA Q 
and GGA [3[ functionals to calculate the CO adsorption 
energies with a numerical uncertainty of ±0.02 eV [lij ]. 
Supercells with relaxed, symmetric five layer slabs are 
employed to model the low-coverage adsorption at f/9 
monolayer on the extended surface. In good agreement 
with previous studies [1, [H, HH , we obtain the fee hol- 
low site to be more stable than the correct top site (by 
0.33 eV within the LDA, by O.lfeV within the GGA). 

Our approach works as follows (we are using here DFT- 
LDA as the starting point, but, of course, one could as 
well start with other functionals, e.g. DFT-GGA): 

1. Do supercell calculations for the extended system 
using DFT-LDA. 

2. Do cluster calculations with the same functional 
and same geometry as in step 1. It may be convenient 
to saturate the cluster surface by hydrogen adatoms, but 
there is in principle no need to do so. 

3. Do corresponding calculations for exactly the same 
cluster as in step 2, but using an improved xc treatment. 
This may employ the B3LYP functional, or Hartree-Fock 
(HF) plus M0ller-Plesset perturbation theory (MP2), a 
coupled-cluster, or a quantum Monte Carlo calculation. 

The difference of the results of steps 3 and 2, i.e. 
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then is the xc correction. Why should this cluster quan- 
tity, iJ xccorr - be of much relevance for the extended sys- 
tem? Figure Q] shows how £i xccorr - changes with cluster 
size using the CO adsorption energy in the fee hollow 
site at Cu(llf ) as example. These calculations [lg, [f7| 
were performed at the DFT-LDA, GGA, and B3LYP 
levels, as well as with HF-MP2. The correction of the 
LDA result is dramatic: It is of the order of several hun- 
dreds eV, and the different xc treatments all give no- 
ticeably different results. In this respect we note that 
none of the employed methods (LDA, GGA, B3LYP, 
HF-MP2) fulfills the variational principle for the ground 
state of the true many-electron Hamiltonian, i.e. they 
all could, in principle, give results that are even be- 
low the true many-body ground state energy. For a 
free CO molecule, for example, the total energies are 
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FIG. 1: Total energy correction £; xccorr - ; se e eq. Q, with re- 
spect to the LDA as function of cluster size and for xc = GGA, 
B3LYP, and HF-MP2 for the adsorption of CO at Cu(lll) in 
the fee hollow site. 
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FIG. 2: The energy correction difference iJ xccorr ' (f cc ) — 
B xccorr '(top). The dash-dotted line at 0.33 eV marks the 
minimum xc correction required to obtain the correct top ad- 
sorption site. The N = oo result for the GGA-LDA correction 
was obtained by supercell calculations (LAPW/APW+lo). 



-3 059.01 eV (LDA), -3 079.83eV (GGA), -3 083.30eV 
(B3LYP), -3 075.89 eV (HF-MP2), and the experimen- 
tal value is -3 084.72 eV. [3] Thus, the order is LDA > 
HF-MP2 > GGA > B3LYP > experiment, were the dif- 
ferences to the experimental values are between — 25.7eV 
(LDA) and -1.4eV (B3LYP). The trend seen in Fig. Q] 
is thus the same as that of the free CO molecule. 

Obviously, £? xccorr - is very different for different xc 
functionals. Interestingly, if we evaluate differences, e.g. 
for different adsorption sites, these differences rapidly ap- 
proach a constant value. Figure [2] shows the difference 
of E xccorr - for CO at Cu(fff) in the fcc and in the top 
adsorption sites. Thus, xc-approximate total energy dif- 
ferences of the extended surface can be corrected through 
higher-level calculations for finite clusters, and the cluster 
size where the xc energy correction term Ai? xccorr - be- 
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TABLE I: Adsorption energies (in eV) for low-coverage CO 
adsorption into different high-symmetry adsorption sites at 
Cu(lll) and Ag(lll). The values at the GGA and B3LYP 
levels are obtained through the xc energy correction scheme, 
using the LDA numbers as reference. The energy of the lowest 
energy structure is taken as energy zero. 



System 


xc 


top 


fee 


hep 


bridge 


Cu(lll) 


GGA 


+0.13 


+0.01 





+0.05 




B3LYP 





+0.21 


+0.26 


+0.20 


Ag(lll) 


GGA 





+0.22 


+0.15 


+0.14 




B3LYP 





+0.42 


+0.43 


+0.35 



comes constant determines the efficiency of the local cor- 
rection approach. For the GGA we also have the supercell 
result at 1/9 monolayer coverage, where the adsorbate- 
adsorbate interaction is negligible, and we added this 
N = oo result as well. Figure [2] demonstrates that 
Ai? xccorr - is converged already for very small clusters (a 
16 Cu atom cluster appears to be sufficient). We em- 
phasize that this applies to the differences and the LDA 
error. For adsorption or reaction energies such clusters 
are by far too small. Apparently the LDA error is even 
shorter ranged than Kohn's nearsightedness concept [13] 
which refers to interaction energies. 

Using the converged values of /\ ^ E xccor ' [ ^ enables us to 
correct the LDA energy of the slab calculations. The 
GGA correction decreases the wrong LDA preference for 
the fee site, but it can not yet change the energetic order 
of the hollow and top adsorption sites. However, at the 
B3LYP level the top site is now more stable by 0.21 ± 
0.03 eV. Interestingly, an almost identical value of 0.28 ± 
0.04 eV is obtained at the HF-MP2 level. This confirms 
the interpretation of earlier B3LYP studies 
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that the main reason for the wrong site preference of 
LDA and GGA functionals is the self-interaction error 
(also present in the free CO molecule). 

TableUgives the energies for all high-symmetry adsorp- 
tion sites at the Cu(lll) and Ag(lll) surfaces, namely 
top, bridge, fee, and hep. For CO/Cu(lll) the top site 
is now the most stable adsorption site, and the optimum 
diffusion energy barrier for the top-bridge-top pathway 
is 0.20 eV. For the other system, CO at Ag(lll), already 
GGA yields the correct top site as the most stable one 
B G3 1 an d it is interesting to verify that a higher-level 
xc treatment does not spoil this description. As shown 
in Table HI the energetic order is indeed not changed at 
the B3LYP level, and the top site is in fact further stabi- 
lized by 0.20 eV, so that the diffusion barrier for the top- 
bridge-top pathway is more than doubled: from 0.14eV 
(GGA) to 0.35 eV (B3LYP). 

The approach is not just applicable to localized per- 
turbations, like an adsorbate or a defect. It can also be 
used to study the bulk cohesive energy. In this case, how- 
ever, it is not possible to express the correction in terms 
of energy differences so that cluster size and edge effects 
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FIG. 3: Cohesive energy E coh - of Cu (cf. eq. ©) for the 
LDA (black line, open squares), and AE coh ' corrections with 
respect to the LDA for GGA, B3LYP, and HF-MP2. The 
N — oo result for the LDA and the GGA-LDA correction was 
obtained by crystal bulk calculations (LAPW/APW+lo). 



cancel. We therefore now write the total energy of the 
many-atom system as sum over contributions assigned 
to the individual atoms, E = ^Zj-Ej. Here E] is the 
energy contribution due to atom I. In the simplest, yet 
physically meaningful approach for metals, Ej is roughly 
proportional to \fc, where c is the local coordination, i.e. 



the number of nearest neighbors (see, e.g. 
our close-packed Cu clusters we then get 



HIH). For 
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where E atom is the calculated total energy of the free 
atom at the given xc level, N the number of atoms in 
the cluster, N c the number of atoms in the cluster that 
are c-fold coordinated, and E coh - corresponds to the co- 
hesive energy expected for an infinite size cluster. Figure 
[3] shows how E coh - changes with cluster size and how it 
depends on the xc functional. 

It can be seen that the correction 
AE coh (GGA — LDA) = S coh (GGA) - £ coh (LDA), as 
well as those for the other xc functionals, is converged 
already for N — 24 clusters with an uncertainty of 
±0.1 eV. Hartree-Fock shows the same convergence 
behavior as B3LYP but at a value of —3.8 eV. For the 
GGA we also performed calculations for the infinite 
crystal, and this gives the difference between the GGA 
and LDA cohesive energies as —1.06 eV, in very close 
agreement to the converged value we get with our xc 
correction scheme. Neef and Doll [22j had obtained 
values of -1.05 eV and -2.04 eV for the (GGA-LDA) 
and the (B3LYP— LDA) correction, respectively, which 
is rather close to our results given in Fig. [3J However, 
this (B3LYP— LDA) correction is too large to match 
experiment. The experimental cohesive energy for Cu 
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is 3.49 eV [19], which differs by -0.83 eV from the LDA 
cohesive energy. Apparently, B3LYP is very bad in its 
description of Cu bulk. Also our cohesive energy at the 
HF-MP2 level (2.82 ±0.1 eV) is significantly smaller 
than the experimental value, but at this point we cannot 
rule out that the HF-MP2 convergence with cluster size 
is different to that of other treatments. In fact, whether 
or not HF-MP2 should work for metals needs a deeper 
theoretical analysis. A more detailed discussion will be 
given elsewhere. [25| 

All results in Fig. [3] were obtained for the LDA lat- 
tice constant. Of course, we could have easily optimized 
the lattice constants for the different treatments (or we 
could have shown the equations of state for the differ- 
ent xc levels). However, this would have complicated 
the graphs due to the intermingling of different effects. 
The focus of the present work is on metal surfaces and 
bulk. However, the methodology proposed in this paper 
for correcting DFT-slab adsorption, cohesive, and diffu- 
sion energies is also applicable to semiconductors (e.g. 
H 2 at Si(001) [26]), or ionic materials (e.g. NaCl). Not 
surprisingly, here the efficiency of the Axe approach is 
even better (for details see [IB)]). The approach had been 
applied by Tuma and Sauer to protonation reactions in 
zeolites (combining DFT and HF-MP2).[27j This work is 
interesting as here the main source of error at the DFT 
level appears to be the lack of the long-range van der 
Waals tails, and by the correction scheme the authors 
could evaluate these van der Waals contributions to ad- 
sorption. 

In summary, we presented a scheme to locally correct 
the total energy (or total energy differences) for errors 
contained in present-day local or semi-local DFT func- 
tionals, e.g. the self-interaction and the lack of van der 
Waals interactions. When looking at appropriate energy 
differences, a smooth and rapid convergence of the xc 
correction with cluster size is observed. This enables us 
to reach convergence at very small cluster sizes or extra- 
polate to N = oo. At these small cluster sizes, that are 
treatable, e.g. with HF-MP2, it is only the xc correction, 
not the total energy, that is converged. The approach 
is demonstrated by computing the energetic order of the 
high-symmetry sites for the low-coverage adsorption of 
CO at Cu(lll) and Ag(lll). The corrections to potential 
energy surfaces obtained with LDA or GGA are found to 
be significant: For CO diffusion at Ag(lll) energy bar- 
riers are changed by more than a factor of two, and for 
CO at Cu(lll) even the topology is altered. Further- 
more, the approach was applied to evaluate the cohesive 
energy of Cu bulk, enabling us to perform HF-MP2 cal- 
culations (via a systematic extrapolation, starting from 
the LDA) for an extended system. In this paper we only 
addressed changes in energy, but for forces the correction 
is straightforward, analogous to eq. |T|). 
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